In general terms, tomography is used for obtaining an image of a cross section of an object in a given plane. X-ray tomography was the first to be developed (in 1970s) and its use is now routine not only in medicine but in some industrial applications as well (internal inspection of mechanical components and flaw detection in materials, for example).
Subsequently, a number of new tomographic methods aimed at industrial processes have emerged, collectively known as process tomography (Williams and Beck, 1995). The aim of these methods, which started to develop in the mid 1980s, is to produce an image of the phase or component distribution in an industrial process using only external sensors and without causing any perturbation to it, as depicted in FIG. 1, which shows a system of process tomography formed by a tank or pipe (A), a data acquisition system (B), and a computer for image reconstruction (C). In other words, process tomography provides a way of ‘looking’ from the outside to inside the process, with no need for physical intrusion or alteration, and this represents a radically different approach to gathering structural information on the process to a global scale, unlike the traditional methods based on local sampling on a certain number of points. The process may occur, for instance, in mixing or stirring vessels, fluidized bed reactors, separator tanks, or a pipeline carrying multiphase flow.
There is a whole range of principles and techniques that can be exploited in process tomography, including electrical methods based on impedance measurement, ultrasound, magnetic resonance, optical methods and those based on ionizing radiation (X- and gamma-rays). Generally speaking, ionizing radiation methods produce images with the highest definition, but are relatively slow. On the other hand, electrical methods yield low-resolution images but are much faster, robust and relatively inexpensive.
In particular with regard to electrical impedance tomography, or electrical tomography for short, there has been a very noticeable progress in the last few years. This type of tomography has two main modalities: capacitance tomography and resistance tomography. In a capacitance tomography system (Beck et al., 1997; Gamio, 1997; Plaskowski et al., 1995), (as depicted in FIG. 2, formed by a sensor (A), a data acquisition system (B) and an image reconstruction computer (C)), normally used with mixtures where the continuous phase is non-conducting, the sensor employed, as depicted in FIGS. 3 and 4, is made of a circular array of electrodes (B) distributed around the cross-section to be examined, and the capacitance between all the different electrode-pair combinations is measured. With the help of a computer and a suitable image reconstruction algorithm, this information is used to create a map (an image) showing the variation of the dielectric constant (or relative permittivity) inside the sensor area, thus providing an indication of the physical distribution of the various components of the mixture. Relative to FIGS. 3 and 4, the electrodes (B) can be located on the outside of a non-conducting pipe (A), in order to simplify sensor construction and avoid direct contact with the process fluids (E). A second external grounded metallic pipe (C) serves as an electric screen and to provide mechanical resistance. The sensor also has two cylindric guard electrodes (D) that are grounded as well.
The British patent GB 2 214 640, issued Sep. 6, 1989, describes an electrical capacitance tomography system that employs a linear back-projection (LBP) algorithm as the image reconstruction method. However, such reconstruction method produces images of a relatively low definition.
Resistance tomography, on the other hand, is aimed at mixtures where the continuous phase is a conductor of electricity (Plaskowski et al., 1995; Williams y Beck, 1995). In this case, the electrodes are installed flush with the inside surface of the pipe (or vessel) wall and in direct contact with the fluids. A number of different excitation current patterns are applied and the resultant voltages are measured. They are then used to construct a map of the conductivity distribution inside the sensor, which reflects the physical distribution of the mixture components.
In principle, electrical capacitance tomography (ECT) has important applications in multiphase flow measurement, particularly gas-oil two-phase flow, which often occurs in many oil wells. The traditional way to quantify the various fluids produced by an oil well is to separate the mixture by gravity in large tanks, prior to measuring each component separately using conventional single-phase flow meters. In the last decades, multiphase flow meters have appeared which allow the quantification of the produced fluids without the need to separate the mixture (Thorn et al., 1997). However, the multiphase meters currently available suffer from an unwanted sensibility to changes in the flow regime, unless they are equipped with flow mixing or conditioning devices that introduce permanent pressure losses (which ultimately translate into energy loss). This limitation should be avoided through ECT as the flow regime could be determined and used to compensate the response of conventional multiphase meters, or, alternatively, it is possible to design a new type of tomographic multiphase meter, based on analyzing series of ECT images from two slightly separated cross-sections of the pipe (Hammer et al., 1997; Plaskowski et al., 1995). Additionally, ECT has potential applications to imaging, monitoring and controlling numerous industrial multiphase processes.
However, so far the main limiting factor to the practical application of ECT has been the lack of fidelity or accuracy of the images obtained using the available image reconstruction methods, giving rise to a need of improved methods as the one introduced in this invention (Yang y Peng, 2003). Simple direct methods like linear back-projection (LBP) yield relatively poor images that only provide a qualitative indication of the component distribution inside the sensor. One said method is described in British patent GB 2214640, issued Sep. 6, 1989.
On the other hand, more sophisticated methods, based on iterative local optimization techniques, generally require one or more regularization parameters whose optimal value depends precisely on the (unknown) image to be reconstructed, apart from the fact that the regularization employed has the effect of smoothing the image contours, making it more diffuse. One said method is described in British patent GB 2329476 A (issued Mar. 24, 1999).
Thus, there is an urgent need of better and more accurate image reconstruction methods. As an example of this, a patent was issued recently covering an image reconstruction method based on the application of artificial neural networks (American patent U.S. Pat. No. 6,577,700 B1, issued Jul. 10, 2003).
The present invention describes new image reconstruction methods based on simulated annealing and genetic algorithms. According to the laws of physics, electrostatics in particular, the ECT sensor (as can be observed in FIG. 4) can be considered as a special case of a system of charged conductors separated by a dielectric medium, the theory of which was first developed by J. C. Maxwell (1873). In our particular case, the sensor electrodes (B) act as the charged conductors while the insulating pipe end (A) the sensor contents acts as the dielectric medium. For an n electrode sensor, the induced electrode charges qi and the electrode potentials vi are related by the following set of linear equations
                                                        q              1                                            =                                                              c                11                            ⁢                              v                1                                                          +                                                              c                12                            ⁢                              v                2                                                          +                                …                                +                                                              c                                  1                  ⁢                                                                          ⁢                  n                                            ⁢                              v                n                                                                                        q              2                                            =                                                              c                21                            ⁢                              v                1                                                          +                                                              c                22                            ⁢                              v                2                                                          +                                …                                +                                                              c                                  2                  ⁢                                                                          ⁢                  n                                            ⁢                              v                n                                                                          ⋮                                                                                          ⋮                                                                                          ⋮                                                                                                                                                                                                              ⋮                                                              q              n                                            =                                                              c                                  n                  ⁢                                                                          ⁢                  1                                            ⁢                              v                1                                                          +                                                              c                                  n                  ⁢                                                                          ⁢                  2                                            ⁢                              v                2                                                          +                                …                                +                                                              c                nn                            ⁢                              v                n                                                                        (        1        )            
where cii are the self-capacitance coefficients (or just self-capacitances for short) of electrode i, while the others, cij, with i≠j, are the mutual capacitance coefficients (or just mutual capacitances) between electrodes i and j. Put in matrix form, equation (1) becomes
                                          [                                                                                q                    1                                                                                                                    q                    2                                                                                                ⋮                                                                                                  q                    n                                                                        ]                    =                                    [                                                                                          c                      11                                                                                                  c                      12                                                                            …                                                                              c                                              1                        ⁢                                                                                                  ⁢                        n                                                                                                                                                        c                      21                                                                                                  c                      22                                                                            …                                                                              c                                              2                        ⁢                                                                                                  ⁢                        n                                                                                                                                  ⋮                                                        ⋮                                                        ⋰                                                        ⋮                                                                                                              c                                              n                        ⁢                                                                                                  ⁢                        1                                                                                                                        c                                              n                        ⁢                                                                                                  ⁢                        2                                                                                                  …                                                                              c                      nn                                                                                  ]                        ⁡                          [                                                                                          v                      1                                                                                                                                  v                      2                                                                                                            ⋮                                                                                                              v                      n                                                                                  ]                                      ⁢                                  ⁢        or                            (        2        )                                q        =        Cv                            (        3        )            
Since the capacitances have the property of reciprocity, i.e., cij≠cji, there are only m=½n(n−1) independent mutual capacitances, corresponding to lower (or upper) triangular matrix of C, and corresponding also to each one of the m different electrode-pairs that can be formed in the sensor.
The value of the mutual capacitances is a complex non-linear function of the conductor system geometry, and of the spatial distribution of the dielectric constant or relative permittivity (hereafter called just ‘permittivity’) of the dielectric medium. In the case of the ECT sensor, the geometry of the electrodes, that of the pipe, and the value of the dielectric constant of the latter, are all fixed. Therefore, it can be said that the mutual capacitances are a function only of the spatial distribution of the dielectric constant inside the sensor, ∈(x,y). The problem of calculating the mutual capacitances corresponding to a specific permittivity distribution inside the sensor is referred to as the forward problem.
The use of the cylindrical end guards (FIG. 3) and the assumption that the phase (and thus the permittivity) distribution does not change too much in the axial direction, allows the sensor to be represented by a two-dimensional (2-D) model (Xie et al., 1989). Unless otherwise stated explicitly, in what follows a 2-D model of the sensor will be used. Therefore, the electric charges qi and the capacitances cii y cij that appear in equations (1) and (2) should be considered quantities per unit length of the electrodes in the axial direction. A tilde (i.e., {tilde over (q)}i, {tilde over (c)}ii and {tilde over (c)}ij) shall be used to denote the total quantities that result from considering the actual length of the electrodes. The previous variables are related between them through the electrode length L, according to
                                          q            i                    =                                                    q                ~                            i                        L                          ,                              c            ii                    =                                                                                          c                    ~                                    ii                                L                            ⁢                                                          ⁢              and              ⁢                                                          ⁢                              c                ij                                      =                                                            c                  ~                                ij                            L                                                          (        4        )            
If the interior of the 2-D sensor is divided into p equal-area regions (or ‘pixels’) where the permittivity is considered constant, then the discrete version of the forward problem is
                    c        =                              [                                                                                c                    1                                                                                                ⋮                                                                                                  c                    m                                                                        ]                    =                                    f              ⁡                              (                ɛ                )                                      =                                                            [                                                                                                                                          f                            1                                                    ⁡                                                      (                            ɛ                            )                                                                                                                                                              ⋮                                                                                                                                                                  f                            m                                                    ⁡                                                      (                            ɛ                            )                                                                                                                                ]                                ⁢                                                                  ⁢                with                ⁢                                                                  ⁢                ɛ                            =                              [                                                                                                    ɛ                        1                                                                                                                        ⋮                                                                                                                          ɛ                        p                                                                                            ]                                                                        (        5        )            
where c is the vector of mutual capacitances (per unit length), ƒi are non-linear functions not known explicitly and ∈ is the vector of permittivities corresponding to the p regions or pixels within the sensing zone.
Applying Gauss's Law, the mutual capacitances per unit axial electrode length can be calculated as
                              c          ij                =                                            q              i                                      v              j                                =                                                    -                                                      ɛ                    0                                    V                                            ⁢                                                ∮                                      Γ                    i                                                  ⁢                                                      (                                          ɛ                      ⁢                                                                                          ⁢                      ▽                      ⁢                                                                                          ⁢                                              ϕ                        j                                                              )                                    ·                                      ⅆ                    l                                                                        =                                          -                                                      ɛ                    0                                    V                                            ⁢                                                ∮                                      Γ                    i                                                  ⁢                                  ɛ                  ⁢                                                            ∂                                              ϕ                        j                                                                                    ∂                      n                                                        ⁢                                      ⅆ                    l                                                                                                          (        6        )            
where ∈o is a physical constant called the permittivity of free space, equal to 8.854×10−12 farads per meter, ∈i is a closed curve surrounding electrode i, dl is a normal vector representing an element of the curve Γi, dl is an element of length of that curve, the symbol ‘•’ represents the scalar product of two vectors, and φj is the electrostatic potential produced in the sensor when applying a voltage of V volts to electrode j (which is called source or excitation electrode) and 0 volts to all others (called detection electrodes).
The potential φj is determined by the solution of the following partial differential equation∇·∈(x,y)∇φj=0  (7)
subject to the boundary conditions (a) φj=V volts on the source electrode and (b) φj=0 on the detection electrodes and the outer screen. In general, equation (7) does not have an analytic solution and must be solved numerically.
The problem of estimating what is the spatial permittivity distribution inside the sensor that corresponds to a specific set of mutual capacitance values is referred to as the inverse problem, and is the problem that image reconstruction methods must address and solve. Normally, the permittivity estimation is made in a discrete way, representing it as a vector ∈ like the one in equation (5), which must be calculated from a vector of observed mutual capacitances c, obtained using of a suitable measurement apparatus.
In order to solve the inverse problem, most ECT systems employ the linear back-projection algorithm (LBP) (Plaskowski et al., 1995; Yang y Peng, 2003; Xie et al., 1989, 1992), which is described next. As a first step, a sensitivity map must be calculated for each one of the m=½n(n−1) possible electrode pairs, given by
                                                        s              i                        ⁡                          (              k              )                                =                                                                                                                c                      i                                        ⁡                                          (                      k                      )                                                        -                                      c                                          i                      ⁡                                              (                        emp                        )                                                                                                                                  c                                          i                      ⁡                                              (                        full                        )                                                                              -                                      c                                          i                      ⁡                                              (                        emp                        )                                                                                                        ⁢                                                          ⁢              for              ⁢                                                          ⁢              i                        =            1                          ,        …        ⁢                                  ,        n                            (        8        )            where k is the pixel number (from 1 to p), ci(k) is the capacitance measured with electrode pair i when the area of pixel k is full of a high-permittivity material while the rest of the sensor is full of a low-permittivity material, whereas ci(full) and ci(emp) are the capacitances for electrode pair i when the sensor is full of high- and low-permittivity material, respectively. Generally, these sensitivity maps are calculated by solving numerically equation (7) and applying equation (6).
Having determined the sensitivity maps, they can be used to obtain a permittivity image from any vector of m=½n(n−1) measured mutual capacitances, c, corresponding to some unknown material distributions inside the sensor. For this, the measured capacitance readings must first be normalized according to
                                          λ            i                    ⁡                      (            k            )                          =                                            c              i                        -                          c                              i                ⁡                                  (                  emp                  )                                                                                        c                              i                ⁡                                  (                  full                  )                                                      -                          c                              i                ⁡                                  (                  emp                  )                                                                                        (        9        )            
where μi is the normalized capacitance for electrode pair i and ci is the actual capacitance measured with that electrode pair.
The basic LBP formula calculates a ‘grey level’ g(k) for each pixel as
                                          g            ⁡                          (              k              )                                =                                                                                          ∑                                          i                      =                      1                                        m                                    ⁢                                                            λ                      i                                        ⁢                                                                  s                        i                                            ⁡                                              (                        k                        )                                                                                                                                  ∑                                          i                      =                      1                                        m                                    ⁢                                                            s                      i                                        ⁡                                          (                      k                      )                                                                                  ⁢                                                          ⁢              for              ⁢                                                          ⁢              k                        =            1                          ,        …        ⁢                                  ,        p                            (        10        )            In principle, this grey level is supposed to be linearly related to the permittivity, with g=1 and g=0 corresponding to the permittivities of the high- and low-permittivity materials, respectively. LBP is based on making a linear approximation to a problem that, as already mentioned, is essentially non-linear (Gamio and Ortiz-Aleman, 2003). Therefore, this image reconstruction method causes considerable errors, which are particularly grave if there are large differences in permittivity in the image.
So far, the main alternative to LBP has been the use of iterative methods that seek to minimize some objective function, employing local optimization techniques like the regularized Newton-Raphson method or other similar approaches (Yang and Peng, 2003). As an example of these methods, there is the one used in the EIT2D software package (Vauhkonen et al., 2001), developed by researchers from Finland and the United Kingdom. Their method is based on minimizing, with respect to ∈, the following functional∥cmeas−ccalc∥2+α2∥L∈∥2  (11)where α is a regularization parameter, L is a ‘regularization’ matrix containing some type of a-priori smoothness information about c, and ccalc=f(∈) is the vector of n calculated mutual capacitance values for a given permittivity distribution inside the sensor. Starting with an initial guess ∈o, the minimization is carried out by the following iterative procedure (basically a Newton-type method with Tikhonov regularization)∈k+1=∈k+[JkTJk+α2LTL]−1{JkT[cmeas−f(∈k)]−α2LTL∈k}  (12)where Jk is the so-called Jacobian matrix of the partial derivatives of f(∈), evaluated at ∈k
                              J          k                =                              J            ⁡                          (                              ɛ                k                            )                                =                                    [                                                                    ∂                                          f                      i                                                                            ∂                                          ɛ                      j                                                                      ⁢                                  ❘                                      ɛ                                          k                      j                                                                                  ]                        ⁢                                                  ⁢                          (                                                i                  =                                      1                    ⁢                                                                                  ⁢                    …                    ⁢                                                                                  ⁢                    m                                                  ,                                  j                  =                                      1                    ⁢                                                                                  ⁢                    …                    ⁢                                                                                  ⁢                    p                                                              )                                                          (        13        )            
However, these image reconstruction methods have the problem that they require, for their correct operation, one or more regularization parameters whose right value is strongly dependent, precisely, on the image that one wishes to reconstruct, implying that one would need to know beforehand the solution to the problem. Moreover, these methods produce distorted images, because the regularization has an excessive smoothing effect on the obtained permittivity. If the regularization is too strong the smoothing effect will occur, and if it is too weak the method can become unstable and/or not converge to the desired solution.
These local optimization algorithms, during their search, explore only a relatively small sector of the solution domain, restricted to the vicinity of the initial guess. If the optimal solution of the problem, i.e., the absolute minimum of the objective function, is located far from the initial guess, it will hardly be reached due to the presence of relative minima in their way, places where these methods can become trapped. The most used methods in this category are least-squares linear inversion and techniques that utilize the gradient of the objective function, like the steepest-descent and the conjugate-gradient methods. In general, local search methods exploit the (scarce) information derived from the comparison of a small number of models (solutions), thus avoiding an extensive search in the whole model space (Sambridge y Drijkoningen, 1992).
Global optimization methods explore the whole solution domain during the inversion process. They carry out an extensive scan within the model space. In this way, despite the existence of partial solutions to the problem, there is a greater possibility that the final solution corresponds to the best fit between the observed data and the synthetic ones. This type of methods, contrary to local techniques, does not require the information provided by the derivatives of the objective function, because in this case the problem does not need to be linearized. Global optimization algorithms use stochastic criteria in order to simultaneously explore all the solution space in search of the optimal model. The best known of the global methods is Monte Carlo, which performs a purely random and unbiased search. In other words, when generating each new model, it does not take advantage of the information obtained from the previously evaluated models (Gallagher et al., 1991). The unguided randomness is the most characteristic feature of this method, which distinguishes it from the rest of the global methods. Among the global optimization techniques, there are also the method of simulated annealing and genetic algorithms. Both were conceived as analogies of optimization systems occurring in nature. Genetic algorithms emulate the mechanisms of biological evolution while simulated annealing is based on thermodynamics. Both methods are inherently non-linear and, therefore, lend themselves naturally to their application in capacitance tomography, a non-linear problem.